##########################################################################################
# 
# Replication Code for:
#
# Carlo Koos & Richard Traunmüller
# "The Gendered Costs of Stigma: 
# How Experiences of Conflict-related Sexual Violence Encourage Civic Engagement
# for Women and Men" American Journal of Political Science
#
# Author: Richard Traunmüller
# Date: 15.08.2023
# Contact: traunmueller@uni-mannnheim.de
#
#
# All analyses were run using:
# R version 4.1.2 (2021-11-01)
# MacOS version 11.4
##########################################################################################

# Clear Work Space
rm(list=ls())

# Load Required Packages
library(list)
library(psych)
library(arm)
library(misreport)
library(stargazer)
library(bayeslist) # Needs to be installed manually! See README.txt for instructions.
library(xtable)
library(denstrip)

# Set Working Directory
# Make sure you set YOURPATH correctly.
#setwd("YOURPATH/AJPS_replication_materials/")
setwd("~/Dropbox/Paper_listexperiment_threecases/Submission/AJPS/AJPS_replication_materials/")

#########################################################################################
# Load and Prepare Data

load("Data/ajps_replication_data_congo.RData") # DRC data

load("Data/ajps_replication_data_liberia.RData") # Liberia data

load("Data/ajps_replication_data_srilanka.RData") # Sri Lanka data

source("Code/ajps_replication_create_variables.R") # Create variables


#########################################################################################
# Prepare List Experiments, Balance Tests, and Tests for No Design Effects

source("Code/ajps_replication_list_experiments.R") # Prepare list experiments


# Table D.3: Balance Test: Covariates regressed on treatment variable.

balance.congo <- summary(lm(treatment ~ female + age.z + edu.z + income.z + hh_size.z + murder_yes + leavehome_yes, data=congo))

balance.liberia <- summary(lm(treatment ~ female + age.z + edu.z + income.z + hh_size.z + cw_kill + cw_displaced, data=liberia))

balance.srilanka <- summary(lm(treatment ~ female + age.z + edu.z + income.z + hh_size.z + killed + displace, data=sri))

xtable(round(cbind(coef(balance.congo)[,1], coef(balance.liberia)[,1], coef(balance.srilanka)[,1]), 2))


# Table E.4: Test of no design effect for three list experiments.

xtable(design.congo$pi.table) # DRC
design.congo$p

xtable(design.liberia$pi.table) # Liberia
design.liberia$p

xtable(design.sri$pi.table) # Sri Lanka
design.sri$p

#########################################################################################
# The Prevalence of Conflict-Related Sexual Violence
# Present Results of List Experiments and Compare to Direct Items

source("Code/ajps_replication_direct_items.R")

# Table G.25 Estimated Prevalence Rates of CSRV
tab <- round(rbind(
  cbind(c.ind.congo, c.ind.liberia, c.ind.sri),
  cbind(se.ind.congo, se.ind.liberia, se.ind.sri),
  
  cbind(c.dir.congo, c.dir.liberia, c.dir.sri),
  cbind(se.dir.congo, se.dir.liberia, se.dir.sri),
  
  cbind(c.diff.congo, c.diff.liberia, c.diff.sri),
  cbind(se.diff.congo, se.diff.liberia, se.diff.sri)), 2)
colnames(tab) <- c("DRC", "Liberia", "Sri Lanka")
rownames(tab) <- c("List Experiment", "", "Direct Questions", "", "Difference", "")
xtable(tab)

# Figure 1: Different estimates of conflict-related sexual violence prevalence with 90% confidence intervals:
# List experiments, direct survey question, and their difference (i.e. degree of sensitivity).

source("Code/ajps_replication_figure_1.R")

#########################################################################################
# G.1 Risk Factors of CRSV and Misreporting

# Table G.5: a) Sensitive item outcome and b) misreporting equations of a multivariate
# regression model of conflict-related sexual violence for DRC.

eady.congo <- listExperiment(Y ~ female + age.z + edu.z + income.z + hh_size.z,
                        data = congo, J = 3,
                        treatment = "treatment", direct = "direct",
                        control.constraint = "none",
                        sensitive.response = 1,
                        misreport.treatment = F)
summary(eady.congo)
cbind(round(eady.congo$par.sensitive, 2), round(eady.congo$se.sensitive, 2))
cbind(round(eady.congo$par.misreport, 2), round(eady.congo$se.misreport, 2))


# Table G.6: a) Sensitive item outcome and b) misreporting equations of a multivariate
# regression model of conflict-related sexual violence for Liberia.

liberia$mon <- ifelse(liberia$direct.1==1 & liberia$Y==0, 0, 1)
eady.liberia <- listExperiment(Y ~ female + age.z + edu.z + income.z + hh_size.z,
                        data = liberia[liberia$mon==1,], J = 3,
                        treatment = "treatment", direct = "direct.1",
                        control.constraint ="full",
                        sensitive.response = 1,
                        misreport.treatment = F)
summary(eady.liberia)
cbind(round(eady.liberia$par.sensitive, 2), round(eady.liberia$se.sensitive, 2))
cbind(round(eady.liberia$par.misreport, 2), round(eady.liberia$se.misreport, 2))


# Table G.7: a) Sensitive item outcome and b) misreporting equations of a multivariate
# regression model of conflict-related sexual violence for Sri Lanka.

sri$mon <- ifelse(sri$direct.1==1 & sri$Y==0, 0, 1)
eady.sri <- listExperiment(Y ~ female + age.z + edu.z + income.z + hh_size.z,
                         data = sri[sri$mon==1,], J = 3,
                         treatment = "treatment", direct = "direct.1",
                         control.constraint = "none",
                         sensitive.response = 1,
                         misreport.treatment = F)
summary(eady.sri)
cbind(round(eady.sri$par.sensitive, 2), round(eady.sri$se.sensitive, 2))
cbind(round(eady.sri$par.misreport, 2), round(eady.sri$se.misreport, 2))


#########################################################################################
# G.2 Risk Factors of Other Forms of Violence

# Table G.8: DRC: Risk Factors of Other Forms of Violence and Direct Item on Conflict-
# Related Sexual Violence.

d.congo <- lm(direct ~ female + age.z + edu.z + income.z + hh_size.z,  data=congo)
summary(d.congo)

vio.1.congo <- lm(murder_yes ~ female + age.z + edu.z + income.z + hh_size.z,  data=congo)
summary(vio.1.congo)

vio.2.congo <- lm(leavehome_yes ~ female + age.z + edu.z + income.z + hh_size.z,  data=congo)
summary(vio.2.congo)

stargazer(d.congo, vio.1.congo, vio.2.congo)


# Table G.9: Liberia: Risk Factors of Other Forms of Violence and Direct Item on Wartime
# Sexual violence.

d.liberia <- lm(direct.1 ~ female + age.z + edu.z + income.z + hh_size.z,  data=liberia)
summary(d.liberia)

vio.1.liberia <- lm(cw_kill ~ female + age.z + edu.z + income.z + hh_size.z,  data=liberia)
summary(vio.1.liberia)

vio.2.liberia <- lm(cw_displaced ~ female + age.z + edu.z + income.z + hh_size.z,  data=liberia)
summary(vio.2.liberia)

stargazer(d.liberia, vio.1.liberia, vio.2.liberia)


# Table G.10: Sri Lanka: Risk Factors of Other Forms of Violence and Direct Item on
# Wartime Sexual violence.

d.sri <- lm(direct.1 ~ female + age.z + edu.z + income.z + hh_size.z, data=sri)
summary(d.sri)

vio.1.sri <- lm(killed ~ female + age.z + edu.z + income.z + hh_size.z, data=sri)
summary(vio.1.sri)

vio.2.sri <- lm(displace ~ female + age.z + edu.z + income.z + hh_size.z, data=sri)
summary(vio.2.sri)

vio.3.sri <- lm(trauma ~ female + age.z + edu.z + income.z + hh_size.z, data=sri)
summary(vio.3.sri)

stargazer(d.sri, vio.1.sri, vio.2.sri, vio.3.sri, type="text")

#########################################################################################
# G.3 Gender Differences in the Civic Effects of CRSV

# Run this code if you want to replicate the MCMC simulations for Liberia
# source("Submission/AJPS/AJPS_replication_materials/ajps_replication_gender_specific_models.R")

# Run this code to use already existing MCMC simulations for Liberia
source("Code/ajps_replication_gender_specific_models_do_not_run_MCMC.R")

# Table G.11: Gender Differences in the Civic Effect of Conflict-Related Sexual Violence:
# DRC

xtable(cbind(round(drc.out.f$par.outcome,2), round(drc.out.f$se.outcome,2), round(drc.out.m$par.outcome,2), round(drc.out.m$se.outcome,2)))

# Table G.12: Gender Differences in the Civic Effect of Conflict-Related Sexual Violence:
# Liberia

xtable(cbind(round(lib.out.f$summaryout[18:26, c(1, 3)], 2), round(lib.out.m$summaryout[18:26, c(1, 3)], 2)))
N.f <- dim(lib.out.f$X)
N.m <-dim(lib.out.m$X)

# Table G.13: Gender Differences in the Civic Effect of Conflict-Related Sexual Violence:
# Sri Lanka

xtable(cbind(round(sri.out.f$par.outcome,2), round(sri.out.f$se.outcome,2), round(sri.out.m$par.outcome,2), round(sri.out.m$se.outcome,2)))

#########################################################################################
# G.4 Results for Civic Engagement w/o CRSV

# Table G.14: DRC

drc.0 <- glm(soc_part.2 ~ female + age.z + edu.z + income.z + hh_size.z + murder_yes 
               + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 , data=congo, family=binomial("logit"))

summary(drc.0)

# Table G.15: Liberia

lib.0 <- glm(outcome_ca ~ female + age.z + edu.z + income.z + hh_size.z + cw_kill 
               + county.1 + county.2, data=liberia, family=binomial("logit"))

summary(lib.0)

# Table G.16: Sri Lanka

sri.0 <- glm(soc_part.2 ~ female + age.z + edu.z +  income.z + hh_size.z +  killed + trauma +
                          + prov.2 + prov.3 + prov.4 + prov.5 + prov.6 + prov.7 + prov.8 + prov.9, data=sri, family=binomial("logit")) 

summary(sri.0)

#########################################################################################
# G.5 Results for Civic Engagement
# Main Results

source("Code/ajps_replication_main_models.R")

# Table G.17: Regression models of civic participation for DRC.

xtable(cbind(round(m.out.congo$par.outcome, 2), round(m.out.congo$se.outcome, 2), 
             round(c(coef(m.out.congo.direct)[-2], coef(m.out.congo.direct)[2]), 2), round(c(se.coef(m.out.congo.direct)[-2], se.coef(m.out.congo.direct)[2]), 2)))

# Table G.18: Regression models of civic participation for Liberia.

xtable(cbind(round(m.out.liberia$par.outcome,2), round(m.out.liberia$se.outcome,2),
             round(c(coef(m.out.liberia.direct)[-2], coef(m.out.liberia.direct)[2]), 2), round(c(se.coef(m.out.liberia.direct)[-2], se.coef(m.out.liberia.direct)[2]),2)))

# Table G.19: Regression models of civic participation for Sri Lanka.

xtable(cbind(round(m.out.sri$par.outcome,2), round(m.out.sri$se.outcome,2),
             round(c(coef(m.out.sri.direct)[-2], coef(m.out.sri.direct)[2]), 2), round(c(se.coef(m.out.sri.direct)[-2], se.coef(m.out.sri.direct)[2]),2)))

# Figure 2: The effect of conflict-related sexual violence on civic participation in three post-conflict
# contexts. Averaged differences in predicted probability along with 90% intervals.

source("Code/ajps_replication_figure_2.R")

#########################################################################################
# Sex Differences in the Civic Effect of Conflict-Related Sexual Violence

# Figure 3: Gender differences in the effect of conflict-related sexual violence on civic engagement in three
# post-conflict contexts. Averaged differences in predicted probability along with 90% intervals.

source("Code/ajps_replication_figure_3.R")

#########################################################################################
# G.6 Results for Donations

# DRC
don.drc <- ictreg.joint(Y ~ female + age.z + edu.z + income.z  + hh_size.z +  murder_yes 
                          + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev,  
                          treat="treatment", 
                          outcome="donate",
                          outcome.reg="logistic",
                          constrained=TRUE,
                          J=3, data=congo) # List Experiment Model
 
don.drc.direct <- glm(donate ~ direct + female + age.z + edu.z + income.z + hh_size.z +  murder_yes
           + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev, data=congo, 
           family=binomial("logit")) # Direct Item Model

# Table G.20: Regression models of donations for DRC
xtable(cbind(round(don.drc$par.outcome,2), round(don.drc$se.outcome,2),
             round(c(coef(don.drc.direct)[-2], coef(don.drc.direct)[2]), 2), round(c(se.coef(don.drc.direct)[-2], se.coef(don.drc.direct)[2]),2)))

# Liberia
don.liberia <- ictreg.joint(Y ~ female + age.z + edu.z + income.z + hh_size.z + cw_kill +  as.factor(county), 
                          treat="treatment", 
                          outcome="trust4",
                          outcome.reg="linear",
                          constrained=TRUE,
                          J=3, data=liberia) # List Experiment Model

don.liberia.direct <- lm(trust4 ~ direct.1 + female + age.z + edu.z + income.z + hh_size.z +  cw_kill 
          + county.1 + county.2, data=liberia) # Direct Item Model

# Table G.21: Regression models of donations for Liberia
xtable(cbind(round(don.liberia$par.outcome,2), round(don.liberia$se.outcome,2),
             round(c(coef(don.liberia.direct)[-2], coef(don.liberia.direct)[2]), 2), round(c(se.coef(don.liberia.direct)[-2], se.coef(don.liberia.direct)[2]),2)))

#########################################################################################
# G.7 Results for Interethnic Relations

source("Code/ajps_replication_interethnic.R")

# Table G.22: Regression models of interethnic relations for DRC.

xtable(cbind(round(int.drc$par.outcome,2), round(int.drc$se.outcome,2),
             round(c(coef(int.drc.direct)[-2], coef(int.drc.direct)[2]), 2), round(c(se.coef(int.drc.direct)[-2], se.coef(int.drc.direct)[2]),2)))


# Table G.23: Regression models of interethnic relations for Liberia.

xtable(cbind(round(int.lib$par.outcome,2), round(int.lib$se.outcome,2),
             round(c(coef(int.lib.direct)[-2], coef(int.lib.direct)[2]), 2), round(c(se.coef(int.lib.direct)[-2], se.coef(int.lib.direct)[2]),2)))


# Table G.24: Regression models of interethnic relations for Sri Lanka.

xtable(cbind(round(int.sri$par.outcome,2), round(int.sri$se.outcome,2),
             round(c(coef(int.sri.direct)[-2], coef(int.sri.direct)[2]), 2), round(c(se.coef(int.sri.direct)[-2], se.coef(int.sri.direct)[2]),2)))


# Figure G.1: The effect of conflict-related sexual violence on inter-ethnic relations in three post-conflict
# contexts. Averaged differences in standard deviations/predicted probability along with 90% intervals.

source("Code/ajps_replication_figure_G1.R")


#########################################################################################
# H Sensitivity Analyses

# H1 No Non-Disclosure
# ATTENTION: Simulations may take a while.
# Figure H.2: Sensitivity to the violation of the ‘No Non-Disclosure’ assumption.

# Run this code if you want to replicate the simulations.
#source("Code/ajps_replication_sensitivity_nondisclosure.R")

# Run this code to use already existing simulations.
source("Code/ajps_replication_sensitivity_nondisclosure_do_not_run_simulations.R")


# H2 Unobserved Confounding
# Figure H.3: Sensitivity to unobserved confounding.

source("Code/ajps_replication_sensitivity_unobserved_confounding.R")


# H3 Sample Selection
# Congo
round(mean(.993 * civic.congo.ind + .067*(-1)), 2)
round(quantile(.993 * civic.congo.ind + .067*(-1), c(.05, .95)), 2)

# Liberia
round(mean(.974 * civic.liberia.ind + .026*(-1)), 2)
round(quantile(.974 * civic.liberia.ind + .026*(-1), c(.05, .95)), 2)

# Sri Lanka
round(mean(.994 * civic.sri.ind + .006*(-1)), 2)
round(quantile(.994 * civic.sri.ind + .006*(-1), c(.05, .95)), 2)

#########################################################################################
# I Alternative Mechanism: Post-traumatic Growth

source("Code/ajps_replication_PTG.R")

# Table I.26: Causal mediation analysis for DRC.
xtable(round(cbind(drc.sv.ptg$par.outcome, drc.sv.ptg$se.outcome), 2)) # Mediator Equation
xtable(round(cbind(drc.ptg.sp$par.outcome, drc.ptg.sp$se.outcome), 2)) # Outcome Equation

# Table I.27: Causal mediation analysis for Sri Lanka.
xtable(round(cbind(sri.sv.ptg$par.outcome, sri.sv.ptg$se.outcome), 2)) # Mediator Equation
xtable(round(cbind(sri.ptg.sp$par.outcome, sri.ptg.sp$se.outcome), 2)) # Outcome Equation

# Figure I.4: Causal mediation analysis: The civic effect of CRSV is not mediated by posttraumatic
# growth (PTG).

source("Code/ajps_replication_figure_I4.R")

##########################################################################################
# END
##########################################################################################




